Toward a systematic 1/d expansion: Two particle properties 
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We present a procedure to calculate 1/d corrections to the two-particle properties around the infinite 
dimensional dynamical mean field limit. Our method is based on a modified version of the scheme 
of Ref. jjj]. To test our method we study the Hubbard model at half filling within the fluctuation 
exchange approximation (FLEX), a selfconsistent generalization of iterative perturbation theory. 
Apart from the inherent unstabilities of FLEX, our method is stable and results in causal solutions. 
We find that 1/d corrections to the local approximation are relatively small in the Hubbard model. 
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During the past few years dynamical mean field the- 
ory (DMFT) became one of the most popular methods 
to study strongly correlated systems 0. DMFT devel- 
oped from the path-breaking observation Q that in the 
limit d — > oo of a d-dimensional lattice model with suit- 
ably rescaled hopping parameters, spatial fluctuations 
are completely suppressed and the self-energy becomes 
local. As a consequence, the self-energy can be writ- 
ten as a functional of the on-site Green's function of 
the electrons and the lattice problem reduces to a quan- 
tum impurity problem, where the impurity is embedded 
in a selfconsistently determined environment. The main 
virtue of this method is that it captures all local time- 
dependent correlations and makes possible to study, e.g. 
the Mott-Hubbard transition or the phase diagram of dif- 
ferent Kondo lattices in detail. 

While in the case of the Mott-Hubbard transition the 
transition seems to be driven by the above-mentioned 
local fluctuations, in many cases correlated hopping Q 
or inter-site interaction effects may play a crucial 
role as well, and while some of these effects can be qual- 
itatively captured by a natural extension of the DMFT, 
others are beyond the scope of it and would only ap- 
pear as 1/d corrections. Furthermore, in order to check 
the quality of the local approximation for a finite dimen- 
sional system of interest, it is very important to compare 
it with the size of the appearing 1/d corrections as well. 

Several attempts have been made to partially restore 
some of the spatial correlations lost in the DMFT. One 
of the most successful ones is the cluster approximation 
proposed by Jarrell et al 0. This method has the ad- 
vantage of being causal, however, it requires considerable 
numerical prowess and it is not systematic in the small 
parameter 1/d. Another method based on the system- 
atic expansion of the generating functional has been sug- 
gested by Schiller and Ingersent 0]. However, despite 
of its technical and conceptual simplicity, this method 
has not been used very extensively because it seemed to 
be somewhat unstable and in some cases gave artificial 
non-causal solutions. 



In the present work we first show, that the method of 
Schiller and Ingersent (SI) can be considerably stabilized 
by a minor, however crucial change in the algorithm, as- 
suring that the contributions of some unwanted spurious 
diagrams exactly cancel. The price for this stability is 
a somewhat increased computation time, since in each 
cycle of the original algorithm an additional subcycle is 
needed to assure cancellation. With this change the SI 
method can then be safely extended to the calculation 
of two-particle properties. Here the main difficulties are 
connected to the inversion involved in the solution of the 
Bethe-Salpeter equation and the non-locality of the irre- 
ducible vertex functions. These difficulties are cicrcum- 
vented by introducing bond variables for the two-particle 
propagators. Finally, we test the general formalism with 
the fluctuation exchange approximation (FLEX) || . 

Although the method presented applies to arbitrary 
lattice structures and various models with nearest neigh- 
bor interactions, for concreteness, let us consider the 
Hubbard model on a <i-dimensional hypercubic lattice at 
half filling: 
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Here the dynamics of the conduction electrons Cj CT is 
driven by the hopping t between nearest neighbor sites, 
rii a = c\ a Ci a is the occupation number, and the electrons 
interact via the on-site Coulomb repulsion U. 

In the SI formalism one considers the following single 
(n = 1) and a two-impurity (n = 2) imaginary time 
effective functionals to generate 1/d corrections: 
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Here, as usually, c a a(j) and c aa (r) denote Grassman 
fields, and the indices a and (3 label the sites for n = 2 
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while they are redundant for n — 1. The 'medium prop- 
agators' QQ-i and must be chosen in such a way that 
the dressed impurity propagators G^ and G^ coincide 
with the full on-site and nearest neighbor lattice propa- 
gators, Ggo* an d Gjf": 

G« = G$ = G 
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In this case one can easily show that — restricting oneself 
to skeleton diagrams of the order of 0(1/ d) — the im- 
purity self energies E^ and E^ 2 2 and the diagonal and 
off-diagonal lattice self energies, E att and E latt are re- 
lated by [| 
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Knowing the lattice self energy the lattice Green function 
can then be expressed as 
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where G® m (z) and G\^(z) denote the unperturbed and 
dressed lattice propagators between sites I and m, respec- 
tively. 

Based on the relations above SI suggested the follow- 
ing simple iterative procedure to obtain a solution that 
includes 0(1/ d) corrections: 

0(1,2) _^ s (l,2) _^ ^latt glatt _^ g(l,2) _ 



A careful analysis shows, however, that the second step 
in this scheme is extremely unstable. To understand this 
it is enough to notice that the second term of Eq. (|J) is 
constructed in such a way that at the fixed point all com- 

(2) 

pletely local skeleton diagrams in the expansion of E-jy 
are canceled by the subtracted E^ term. However, this 
cancellation only happens under the condition that the 
dressed Green's functions G^ and G^ are exactly the 
same. If G$ and G (1) differ by a term of 0(1/ d) in a 
given iteration step, the cancellation above is not exact, 
and an error of the order of 2d x 0(1/ d) ~ 1 is gen- 
erated immediately. Moreover, the generated erroneous 
term is typically acausal because of the subtraction pro- 
cedure involved in Eq. (Q), and may drive the iteration 
towards some more stable but unphysical fixed point of 
the integral equations. We suggest to replace the critical 
steps C/t 1 - 2 ) E (1 ' 2 ) =$> E latt by the following procedure: 
(1) Calculate G^ from £( 2 ) , (2) Determine selfon- 

sistently in such a way that G^ = G^ be satisfied, (3) 
Determine E^ 1 ' 2 ) and from them E latt . Step (2) above 
is crucial to guarantee that unwanted terms in Eq. (^) 
exactly cancel. 

The two-particle properties can be investigated in a 
way similar to Ref. [Q . To this end we introduce the lat- 
tice particle-hole irreducible vertex function r latt . which 
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FIG. 1. Graphical representation of the Bethe-Salpeter 
equation: Double lines indicate dressed Fermion propagators, 
while boxes denote particle-hole irreducible vertex diagrams. 



is connected to the full lattice propagator L latt by the 
Bethe-Salpeter equation (see Fig. pi): 



L latt (w) = V Q att (uj)(l + T latt (uj)V att (oj)) 
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where ui denotes the transverse frequency and a ten- 
sor notation has been introduced in the spatial, spin 
and other frequency indices: H^)ilf 2 % 3 ; Q A , U3 -> 
L(cj). The 'vertex-free' propagator L att (a;) is defined 

as ^{f tt ( Ci; )i 1 Vi2^i;i3, 4 i4,W3 = ^lff3^2^4^iw 3 ^1^3(^1) 

G' a « (i(wi + w)). 

A detailed analysis shows that up to 1/d order the 
only non-zero matrix elements of r latt (ci;) are those where 
the indices i±, 12, 13, ii belong to the same or two nearest 
neighbor lattice sites, i.e. a bond. A thorough investiga- 
tion of the corresponding skeleton diagrams shows that 
r latt (o;) can be expressed similarly to Eqs. (^) and (||) 
as: 
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where in the second line it is implicitely assumed that the 
ifc's are not all equal. Here the particle- hole irreducible 
one- and two-impurity vertex functions, T^ 1 ' , and r*- 2 -* 
are defined similarly to r latt , and satisfy the impurity 
Bethe-Salpeter equation: 



LW( W )=LW( W )(l+fW( W )L( w) ( W ))- 
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with n — 1,2. Of course, in the impurity case the spa- 
tial indices of the propagators Lg™^(w) and L^(w) are 
restricted to the impurity sites, but apart from this the 
Lq (w)'s are equal to the lattice propagator L att (tj). 

From the considerations above it immediately follows 
that the 1/d corrections to the two-particle properties can 
be calculated in the following way: (1) Find the solution 
of the single particle iteration scheme, (2) Determine the 
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FIG. 2. The generating functional for FLEX: Double lines 
denote single- and two-impurity dressed cavity propagators 
G (1) (t,t ) and G (2) (t,t'). Position and spin indices are not 
shown. The corresponding self-energy diagrams are obtained 
as S(t — t') = —8<P[G(t,t')]/8G(t',t), and the particle-hole 
irreducible vertex is given by a similar second order functional 
derivative. 
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one- and two-impurity correlators, (3) Invert Eq. (^|) to 
obtain f W and f ( 2 \ (4) Calculate f latt using Eq. (§), 
and (5) Solve the Bethe-Salpeter equation (^) for L 
and calculate two-particle response functions from it. 

The major difficulties in the procedure above are as- 
sociated with the inversion appearing in Eqs. (Q), since 
the propagator L{, att connects any four lattice sites and 
has an infinite number of frequency indices. The first dif- 
ficulty can be resolved by observing that r latt connects 
only neighboring sites. Therefore with a introduction 
of bond variables a partial Fourier transformation can be 
carried out in these, and summations over all pairs of lat- 
tice sites reduce to a summation over d 'bond direction 
indices' and two additional indices specifying the posi- 
tion of the electron and the hole within a given bond. 
Furthermore, to avoid overcounting, the vertex function 
piatt mug ^ replaced by a slightly modified 'bond vertex 
function' || . A further reduction of the matrices involved 
can be achieved by diagonalizing the propagators in the 
spin labels. Finally, to carry out the summations and in- 
versions over the infinite omega variables we introduced 
a frequency cutoff lu c and extrapolated the lo c = oo re- 
sult from a finite size scaling analysis in this cutoff [^0) , 
thereby reducing the numerical error of our calculations 
below one percent. 

To test the procedure above one is tempted to try to 
generalize the iterative perturbation theory (IPT) ap- 
plied remarkably successfully for the d = oo case [ 111] , 
however, it is clear from the discussion above that within 
IPT it is impossible to satisfy the condition G^ = 
(which explains why earlier attempts to generalize IPT to 
order 1/d failed p2|). We therefore applied the so-called 
fluctuation exchange approximation (FLEX) ||. While 
this method is unable to capture the metal insulator tran- 
sition, it is able to reproduce the Kondo resonance in the 
metallic phase fl3}| , has been successfully used to cal- 
culate weak and intermediate coupling properties of the 
2-dimensional Hubbard model [Of , and it has the impor- 



FIG. 3. Top: Calculated diagonal and off-diagonal lat- 
tice self-energies for (7 = 4 and T — 0.002. All energies are 
measured in units of t. Bottom: Local spectral functions for 
U = and U = 4 at T = 0.002. 

tant property of being formulated in terms of the dressed 
single particle Green's functions. In this approach the in- 
teractions between particles are mediated by fluctuations 
in the particle-particle and particle- hole channels, and 
the self-energies and the particle-hole (particle-particle) 
irreducible vertex functions are generated from the gener- 
ating $ functionals built in terms of the dressed Green's 
functions, depicted in Fig. ||. A further advantage of 
FLEX is that due to the special structure of the dia- 
grams involved a fast Fourier transform algorithm can 
be exploited to increase the speed and precision of the 
calculation substantially. 

The calculated three-dimensional diagonal and off- 
diagonal lattice self-energies are shown in Fig. |^ together 
with the local spectral function. These have been ob- 
tained by means of a Pade approximation to carry out 
the analytic continuation from the imaginary to the real 
axis. Though in the spectral function a well-developed 
Kondo peak is observed, the FLEX is unable to repro- 
duce the depletion of spectral weight in the neighbor- 
hood of it due to the 'over-regularization' characteristic 
to most selfconsistent perturbative schemes. Remark- 
ably, we experienced no convergence problems similar to 
those of Ref. H , apart from the ones inherent in FLEX 
S. We checked that the spectral functions integrate to 
one within numerical precision and the solutions obtained 
are causal. The typical values of are nearly an or- 
der of magnitude smaller |l4| than Sff*', indicating that 
the local approximation gives surprisingly good results 
and 1/d corrections are indeed small as anticipated in 
Ref. Q and also in agreement with the results of Ref. [Q . 
To get some further information about the quality of lo- 
cal approximation in Fig. we plotted the momentum 
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FIG. 4. Momentum dependent spectral functions along 
the line k — a(n,n,n) within the local (dashed lines) and 
1/d calculations (continuous lines). The alpha values used 
were a = 0.497, a = 0.495, and a = 0.49 from left to right. 
The energy of the quasiparticles has been renormalized ap- 
proximately by a factor of four compared to the bare electron 
energies. 



dependent spectral functions at different points of the 
Brillouin zone. The 1/d contributions give typically a 
10-20 percent correction, but none of the generic proper- 
ties is modified in the paramagnetic phase. 

Once convergence is reached at the single particle level, 
one can turn to the two-particle properties. Within 
FLEX this is somewhat simpler, because — although 
many rather complicated diagrams are generated [|S|]9[] 
— r( n ) can be built up directly in terms of the full lat- 
tice Green's functions. We find that similarly to the off- 
diagonal self-energy the off-diagonal elements of r latt are 
rather small. Having solved Eq. (Q) one can calculate 
various correlation functions. As an example, in Fig. |E| 
we show the momentum dependent susceptibility of the 
half-filled Hubbard model in its paramagnetic phase for 
two different temperatures along the (1,1,1) direction, 
obtained from the FLEX calculations. The susceptibility 
develops a peak at (n, tt, tt) at low temperatures, as a sign 
of unstability toward antiferromagnetic phase transition. 

We also determined the transition temperature at sev- 
eral values of U and compared our results with existing 
Monte Carlo data ]l5|] . We found a critical tempera- 
ture T c typically by a factor of three lower than that of 
Rcf . J^] . This difference is a result of the overregulariza- 
tion of the interaction vertex by FLEX. Indeed, replacing 
r latt (w) by the bare particle-hole vertex in Eq. (j7j) the 
order-parameter fluctuations become larger (see Fig. ||) 
and T c is in excellent agreement with the Monte Carlo 
data. 

In conclusion, we presented an extended version of the 
SI method to calculate 1 /d corrections to the two-particle 
properties. We tested the new procedure by FLEX. No 
convergence problems and no violation of causality ap- 
peared in our method, although this is not generally guar- 



FIG. 5. Momentum dependent susceptibilities. The sus- 
ceptibility develops a peak at (tv, tv, tt) as a precursor of the 
antiferromagnetic phase transition. 



anteed within the present scheme. Our method should 
be tested on other models and with other, more time- 
consuming methods in the future as well. 
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